clc
clear

syms lambda theta alpha beta eta R Ra Rb Rr omega idom

panel = [1,5,6];

% for iPanel_type = 1:3
iPanel_type = 2;
    iPanel = panel(iPanel_type);
    if iPanel < 5
        lambda = alpha + ( idom - 1 ) * pi / 2;
        theta  = atan2( cos( alpha ) * tan( beta ), 1 );
    elseif iPanel == 5
        lambda = atan2( -tan( alpha ) / tan( beta ), 1 );
        theta  = atan2( 1 / sqrt( tan(alpha)^2 + tan(beta)^2 ), 1 );
    elseif iPanel == 6
        lambda = atan2( tan( alpha ) / tan( beta ), 1 );
        theta  =-atan2( 1 / sqrt( tan(alpha)^2 + tan(beta)^2 ), 1 );
    end
    
%     R = eta;
%     R = str2sym('R(alpha,beta,eta)');
%     R = str2sym('R(eta)');
    
%     v = [alpha;beta;eta];
    v = [alpha;beta;R];
    
    a11 = R * cos(theta) * diff(lambda,v(1));
    a12 = R * cos(theta) * diff(lambda,v(2));
    a13 = R * cos(theta) * diff(lambda,v(3));
    a21 = R * diff(theta,v(1));
    a22 = R * diff(theta,v(2));
    a23 = R * diff(theta,v(3));
%     a31 = diff(R,v(1));
%     a32 = diff(R,v(2));
%     a33 = diff(R,v(3));
    a31 = 0;
    a32 = 0;
    a33 = 1;
    
    A = [a11,a12,a13;...
         a21,a22,a23;...
         a31,a32,a33];
    iA = inv(A);
    G  = transpose(A) * A;
    iG = inv(G);
    
    disp('Calculate Christoffel Symbol')
    CS = Christoffel_Symbol(3,G,iG,v);
    
    OmegaSph  = omega * [0;cos(theta);sin(theta)]; % The angular velocity of the Earth in spherical coordinate
    OmegaCube = iA * OmegaSph; % The angular velocity of the Earth in cube coordinate
    
    omega2 = omega^2;
    CentrifugalForceSph  = omega2 * [0;-omega2*sin(theta);omega2*cos(theta)]; % Centrifugal Force of the Earth in spherical coordinate
    CentrifugalForceCube = iA * CentrifugalForceSph; % Centrifugal Force of the Earth in cube coordinate
    
    A                    = simplify(A);
    iA                   = simplify(iA);
    G                    = simplify(G);
    iG                   = simplify(iG);
    CS                   = simplify(CS);
    OmegaSph             = simplify(OmegaSph);
    OmegaCube            = simplify(OmegaCube);
    CentrifugalForceSph  = simplify(CentrifugalForceSph);
    CentrifugalForceCube = simplify(CentrifugalForceCube);
    
    fid = fopen('CS.txt','wt');
    for i = 1:3
        for j = 1:3
            for k = 1:3
                fprintf(fid,'%s %s %s \n',str2sym(['CS(',num2str(j),',',num2str(k),',',num2str(i),')']),'=',CS(j,k,i));
            end
        end
    end
    fclose(fid);
    
%     disp('Write A matrix in c language code')
%     ccode(A,'file',['Cube_A_P',num2str(iPanel),'.c']);
%     disp('Write Christoffel Symbol in c language code')
%     ccode(CS,'file',['Cube_Christoffel_Symbol_P',num2str(iPanel),'.c']);
%     ccode(OmegaCube,'file','Cube_Omega.c');
% end

function CS = Christoffel_Symbol(n,G,iG,v)
% Calculate Christoffel symbol
% n dimension
% G metric tensor
% iG inverse metric tensor
% v coordinate vector
CS = sym(zeros(n,n,n)); % (Gamma_jk)^i
for m = 1:n
    for k = 1:n
        for j = 1:n
            for i = 1:n
                CS(j,k,i) = CS(j,k,i) + 0.5 * iG(i,m) * ( diff(G(j,m), v(k)) ...
                                                        + diff(G(k,m), v(j)) ...
                                                        - diff(G(j,k), v(m)) );
            end
        end
    end
end
end



